Automatic Differentiation of C++ Codes for Large-Scale Scientific Computing

We discuss computing first derivatives for models based on elements, such as large-scale finite-element PDE discretizations, implemented in the C++ programming language. We use a hybrid technique of automatic differentiation (AD) and manual assembly, with local elementlevel derivatives computed via AD and manually summed into the global derivative. C++ templating and operator overloading work well for both forwardand reverse-mode derivative computations. We found that AD derivative computations compared favorably in time to finite differencing for a scalable finite-element discretization of a convection-diffusion problem in two dimensions. Computing derivatives is ubiquitous in scientific computing; examples include algorithms for nonlinear equation solving, optimization, stability analysis, and implicit time integration. Computing derivatives quickly and accurately improves both the efficiency and robustness of these numerical algorithms, particularly in the presence of ill-conditioning. In this paper, we discuss computing first derivatives of element-based models implemented in ANSI/ISO C++. We use the term “element” in a broad sense to encompass any model whose computation consists of repeated evaluations of a small set of functions, each involving relatively few of the variables of the overall problem. Many classes of models fall into this category, including finite-element and finite-volume PDE discretizations and network models. We use a hybrid technique of automatic differentiation (AD) and manual assembly similar to [1, 2] to carry out the model evaluation and derivative computation one element at a time. This decomposition is discussed in more detail in Section 1, which generalizes the ideas in [2] to general element-based models and additionally describes how to compute the global adjoint. We focus on ANSI/ISO C++ codes because much modern scientific code development is done in C++. Since no source transformation tools for C++ were available to us, we used C++ operator overloading to implement AD for computing the element-level derivatives. We assume the reader is familiar with AD and the methods for implementing it; see [3] for a good introduction to these concepts. We used two separate AD packages: the public domain package Fad [4] for forward-mode AD and our own reverse-mode package Rad [5]. ?? Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy under Contract DEAC04-94AL85000. This is SAND2006-0902C.